Transient structures in a granular gas 
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A force-free granular gas is considered witli impact-velocity dependent coefficient of restitution 
as it follows from the model of viscoelastic particles. We analyze structure formation in this system 
by means of three independent methods: Molecular Dynamics, numerical solution of the hydrody- 
namic equations and linear stability analysis of these equations. AH these approaches indicate that 
structure formation occurs in force-free granular gases only as a transient process. 
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Cluster and vortex formation in force-free granular 
gases is a most striking phenomenon which makes these 
systems so distinct from gases of elastic particles, like 
molecular gases. First detected and explained by Gold- 
hirsch and Zanetti ^J] and McNamara 2], clustering and 
later vortex formation have been intensively studied 
(e.g. 0,0,0)- Clustering has been also reported in re- 
cent simulations of the hydrodynamic equations 6]. In 
all these studies the simplifying assumption of a constant 
coefficient of restitution e has been used. 

It is known, however, that the coefficient of restitution 
is not a material constant, but a function of the impact 
velocity, which reads for the most simple model of vis- 
coelastic spheres 0: 

e(5)-l-7l5l'/' + (3/5)7'|5l'^'T--- (1) 

Here g is the normal component of the impact velocity, 
and 7 is a known function of the particles' material prop- 
erties 0]. The typical value of e, corresponding to the 
thermal velocity vt — \/'^T{t)/m {T{t) is the tempera- 
ture of the gas and rn = 1 is the particle mass) behaves as 
\ — '-yT^/^'^, i.e., it tends to the elastic limit, e ^ 1, 
as the gas cools down, T — > 0. 

Since a force-free gas of elastic particles tends to be ho- 
mogeneous, one can naively assume that a granular gas of 
viscoelastic particles tends to be finally uniform as well. 
However, this is not necessarily true: The collisions be- 
come perfectly elastic only in the limit 5 = when all par- 
ticles are at rest. If the gas cools down too fast the resid- 
ual structures may get frozen and persist due to a lack 
of kinetic energy. We illustrate this for a more general 
model of e{g), with (1 — e^) ^ . The cooling rate of 
such gas is estimated as T -n(l-e|,)r^/2 ^ -nT^/^+'^ 
P,0, where n is the gas number density. The gas density 
decreases with time due to cluster growth. If we assume 
n ~ t~°' (0 < a < 1) 10, we obtain the estimate for the 
gas temperature, T ^ t~^, z = 2(1 — a)/{l + 2/3), and for 
the gas pressure, P = nT ^ t^"^^ . If Tc\ and rid are re- 
spectively the temperature and number density inside the 
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clusters, similar estimates yield: Tci ^ ~"ciTci ■ 
suming that rid keeps approximately constant, we obtain 



Tc\ ^ t 2/(1+2/3) estimate the pressure in clusters, 
-Pel ^ .f{nci)Tci ^ t-2/(i+2/3)^ where the factor / does not 
depend on T. The necessary condition for clusters to dis- 
solve reads Pel > P for t ^ oo, i.e., a + z > 2/{l + 2/3). 
With z given above this is equivalent to the condition 
(3 > 1/2, that is, only gases with a pronounced depen- 
dence of Et on T evolve to a uniform final state. For 
gases of viscoelastic particles (3 = 1/10 < 1/2, hence it 
is not clear whether structures arise only temporarily in 
these systems, or appear and get frozen. 

To address this problem we study the evolution of a 
gas of viscoelastic particles by means of three indepen- 
dent methods: event-driven Molecular Dynamics (MD), 
numerical solution of the hydrodynamic equations and 
linear stability analysis of the HD equations. 

The MD simulations with periodic boundary condi- 
tions were performed for a 2d gas of iV = 10^ particles 
of diameter a = 1, with the coefficient of restitution ac- 
cording to Eq. The gas has the packing fraction cj) — 
mra'^ /A = 0.2, initial temperature Tq = 1 and the dissipa- 
tive coefficient 7 = 0.0577 which corresponds to the ini- 
tial coefficient of restitution £7- ~ 1 — 7 ~ 0.94. Starting 
the simulations with homogeneous distribution of parti- 
cles (homogeneous cooling state, HCS), clusters appear 
and grow until they reach the system size, then the clus- 
ters dissolve (Fig. Figure |21 shows the evolution of the 
average particle energy, E{t) — {^/'2N)^-mvf, which 
evolves in the HCS according to the modified Haff law 

Eit)/Eo = Th(t)/To = (1 + t/To)-'^^ , (2) 

where 

To-i = {24/5)qoT-\0)5, S = (7/Ci)(ro/m)i/i" . (3) 

Here T^^{t) — 271^92(17) y7rT(t)/TO is the mean collision 
time, 52(0") = (1 ^ 7(j)/16)/{l — 0)^ is the contact value of 
the pair correlation function. The constants qo ~ 0.173 
and Ci « 1.1534 are also known analytically At 
later times the deviation from the Haff law becomes pro- 
nounced, however, as the gas further evolves the clus- 




FIG. 1: MD simulation of A'' = 10^ viscoelastic particles. 
Clusters appear as transient structures. The snapshots where 
taken after 0, 200, 800, 7,500, 20,000 and 70,000 collisions per 
particle. 




FIG. 2; The energy decay as obtained by MD (left) and by 
the numerical solution of the HD equations (right). The dash- 
dotted lines show the modified Haff law, Eq. JJJ. The end of 
the HCS is indicated by the vertical dashed lines. 



ters and vortices dissolve, and the system approaches the 
regime of the homogeneous cooling, see Fig. 

To explain the observed effects we consider the hy- 
drodynamic equations 0, 0| for the number density 
n(r, t), flow velocity u(r, t) and temperature field T(r, t) 

dn ^ . . 

— = -V, (nu,) (4) 
du 1 

-rr^ = - (ujVj) Ui H Vj {riijki^kUi - PSij) (5) 

dT PI 

- (ujVj) T (VjUi) + -r]ijki (VfcU/) (S/jUi) 



dt 



n n 



(6) 



Here P = nT[l + (1 + £)0g2(o-)] and r^ija = tli^ikSji + 
Sii5jk — 5ij5ki) is the viscosity tensor, rj and n respectively 
are the coefficients of viscosity and thermal conductivity; 
C is the cooling coefficient and the coefficient ji is spe- 
cific for granular gases The coefficients rj, k, ^ 
and C. have been recently derived for a gas of viscoelas- 
tic particles |lOl |. They may be written as an expansion 
(6 = ?7, K, fi, C), 




biS' 



(7) 



FIG. 3: The gas density obtained from the numerical solution 
of the HD equations 14I5I6II with periodic boundary and ran- 
dom initial conditions. The parameters 7, 0, To are the same 
as in the MD simulations. Density inhomogeneities (clusters) 
appear but eventually dissolve. The snapshots were taken af- 
ter 0, 200, 800, 7,500, 20,000 and 70,000 collisions per particle, 
in correspondece to the snapshots in Fig. The number of 
collisions was computed using the average temperature. 



where 5'{t) = S [2T{t) /To]'^^^° is the time-dependent dis- 
sipation parameter and the microscopic expressions for 
bo, bi, 62 are given in We wish to stress that the 

temperature dependence of these coefficients differs dras- 
tically from the case e = const, see e.g. \l2, 13, f4j. 

We numerically solve the hydrodynamic equations 
(|4I5I6|I with the coefficients Eq. ^ calculated for the 
same microscopic parameters 7, cr, m as used for the MD 
simulations. We use a finite volume discretization scheme 
of global second order on a staggered grid. The integra- 
tion in time is done through a TVD multi-step scheme of 
fourth order . We use a 50 x 50 computational domain 
with periodic boundary conditions. Special care has been 
taken treating the advection terms in (|4I5I6|I 0| . 

We start the numerical integration with random initial 
conditions for the density and flow- velocity field (thermal 
fluctuations) and conflrm the transient character of the 
pattern formation, Fig. |21 The cooling curve for E{t) 0| 
demonstrates qualitatively the same behavior as observed 
in the MD simulations, Fig. [5] (right). 

To study the mechanism of pattern evolution more 
directly we considered initial conditions with a super- 
imposed sinusoidal mode and observed similar transient 
structures. For reasons explained below, we have chosen 
the transverse velocity (shear) mode. It has the compo- 
nents, Uy{f,0) = Uky{0) sm{k x), Ux{r,0) = 0, T(r, 0) — 
To, n{f, 0) = Uh, where k = {li^jV)!, (/ = 1, 2 . . .) is the 
wave-number, L is the system size and Mfej,(0) is the ini- 
tial amplitude of the mode. Solving the hydrodynamic 
equations numerically we analyze the evolution of the 
shear mode, Fig. ^ the longitudinal mode and the den- 
sity mode, Fig. [S] 
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FIG. 4: Evolution of the reduced amplitude of the shear mode, 
Wfc,±(r) = Uky{T) /vt{t). Initially only a single shear mode 
with k — 2tt/L and with an amplitude of the order of the 
thermal fluctuations is excited. Full lines - numerical solution 
of the hydrodynamic equations, dashed lines - results of the 
linear analysis, Eq. 1114^ . The parameters are the same as 
in Fig. 13 Dimensionless time and length are used (see text). 
Note that the different k correspond to a different system size. 




FIG. 5: Evolution of the reduced amplitude of the longitudi- 
nal mode, u)feji(r) = Ukx{T)/vT{r) (left) and of the density 
mode Pkir) (right). Parameters, units and initial conditions 
are the same as in Fig. 0] 



As it follows from Fig. 0] the shear mode (in reduced 
units) is unstable for wave- vectors smaller than a thresh- 
old value. This instability leads to formation of vortices 
. The growing shear mode initiates the growth of the 
mode of longitudinal velocity which, in its turn, causes 
a growth of the density mode corresponding to cluster- 
ing, see Fig. [S] This indicates that clustering is driven 
mainly by a nonlinear mechanism 0,0: As explained be- 
low, a single transverse mode, cannot excite by a linear 
mechanism, neither a longitudinal mode, nor a density 
mode. Hence only nonlinear coupling between the shear 
and longitudinal modes may cause the observed excita- 
tion of the longitudinal mode |l8l | , which further initiates 
the growth of the density mode, i.e. clustering. Similar 
mechanisms govern the initial stage of pattern formation 
in a gas with e = const P, 0. Contrary to the case 
e = const, in a force- free gas of viscoelastic particles all 
modes eventually decay, and only transient structures ap- 
pear. 

To obtain quantitative estimates we perform a stabil- 



ity analysis of the HD equations H4I5I6(I with respect to 
the HCS at density rih, and temperature Th{t). We as- 
sume that the deviations 5T{r,t) — T[r,t) — Ih(i) and 
5n{f, t) — n(r, t) — rtjj are small, 



ST 



< 1. 



Sn 



< 1 



w 



Vt 



«1 (8) 



and linearize these equations. Below we will use the di- 
mensionless time r, measured in units of Tc(t)/2, and 
dimensionless length f, measured in units of 1^/2 {Iq^ — 
V2nnag2{cr) is the mean free path) and write the lin- 
earized equations for the Fourier transforms of the fields 
9{r,t), p{f) and w{r,t): 



dr 
dr 
dr 



dr 



(9) 
(10) 



[Ak^ + {QI3/^)5'{T)\ei^~ikw^ (11) 

(12) 



Here and wj:^ are respectively the longitudinal (i.e. 

parallel to the wave-vector k) and transversal (i.e. per- 
pendicular to k) components of the Fourier mode w^{t), 
j3 = w 0.323 and fii = 1.811 are the constants. 

To obtain Eqs. ©-(CSl we keep only leading-order terms 
with respect to the dissipative parameter 5'{t) in the ex- 
pansion Eq. (|7|) for the coefficients 77, k, jjL, C, 0. Note 
the important difference between the linearized equations 
(P|)- (tT^ and the corresponding equations for the case 
e = const: In the latter case the coefficients in these 
equations are constant (e.g. |3jI3)j while in the former 
case they depend on time via 



5\t) = 2i/i",5[l-f2go<^T/5]" 
The solution for wj:^{t) reads 



(13) 



(14) 



where wj:^{Q) is the initial amplitude of the shear mode. 
There exists a critical wave-number, k*^ = ^/2qo6, which 
separates two regimes: Shear modes with k> k*^ always 
decay, while those with k < k*^ initially grow and reach 
a maximum, 



" 290(5 ■ 
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^ 5fc2 \ 


fc^ e 


exp 1 





(15) 



at T^f^ ~ 5/fc^ — 5/(2go<5), then they decay and die off 
completely. The formation of vortices is attributed to 
the growth of the shear mode 0,0, therefore, the 
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vortices of size ~ decay after a transient time ~ 
. Any system of size ~ L has a minimal wave-number 
~ 271 /L, hence all shear modes wih decay by the time 
T* - 5LV47r2 - 5/(2(jo(5)- 

To perform the stabihty analysis for the other three 
modes we write Eqs. (I1UII - H12II in the form 



(16) 



The matrix Mfc(T) has time- dependent eigenvalues and 
eigenmodes, which for small dissipation 5 are analogous 
to the sound and the heat mode of a gas with a con- 
stant e 0, • For large wave- numbers k all the modes 
decay, while for small k the heat mode may grow. The 
critical k may be found from the condition = 0, or 
det|Mfc(T)| =0: 



4^/5 



(2. 



/lO 



(35 



1/2 



1 



-1/2 



(17) 



Modes with k > k^^^ (0) always decay, while those with 
fc < (0) may initially grow. Since fcj^ (r) decreases with 
time even a mode with fc < /cj^ (0) which initially grows, 
starts to decay after a transient time, when the condition 
k > fc[^(r) is fulfilled. 

The value of k*^ (r) becomes smaller than the minimal 

wave-number 2tt/L at time T|* ~ — 5/{2qo6), i.e., 

for T > T|* the amplitude of any of the modes pj, 
decays. Note that Jj\ agreement with the sim- 
ulation results. Figs. 01 [sl Since clustering and vortex 
formation is attributed to the instability of the heat and 
shear modes P, 0, 0| we conclude that the eventual de- 
cay of all modes predicted by the linear stability analysis, 
implies the transient structure formation. 

We have studied a force-free granular gas of viscoelas- 
tic particles by means of MD, numerical solution of the 
HD equations and linear stability analysis of the HD 
equations. All three methods indicate that structure for- 
mation in a gas of viscoclastic particles occurs only as 
a transient phenomenon, whose duration increases with 
the system size N for the same particle number density. 
Correspondingly, for larger N the structures appear to 
be denser and more long-lived, whereas an extrapola- 
tion to the unbounded system cannot be easily concluded 
from these arguments. The finite duration of the cluster 
phase is in a sharp contrast to the case of a gas with 
a simplified collision model e = const, where structures 



have been proven to arise and to continuously develop. 
In our simulations, due to the limited system size and 
periodic boundary conditions, we have unphysical self- 
interactions of the clusters. We believe, however, that 
this effect does not invalidate the main conclusion of our 
study of the transient character of structure formation in 
force-free granular gases. 
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